

#---Portmann & Stojanovic 2017: Electoral discrimination against immigrant-origin candidates---#

#load("PrefVotes")

#-------------------- ESTIMATE MODEL

d$CodeBerufK32 <- ifelse(d$CodeBerufK3 == 2, 1, 0)
d$CodeBerufK33 <- ifelse(d$CodeBerufK3 == 3, 1, 0)
d$party_c32 <- ifelse(d$party_c3 == 1, 1, 0)
d$party_c33 <- ifelse(d$party_c3 == 2, 1, 0)
d$communeAdliswil <- ifelse(d$commune == "Adliswil", 1, 0)
d$communeBuelach <- ifelse(d$commune == "Buelach", 1, 0)
d$communeDietikon <- ifelse(d$commune == "Dietikon", 1, 0)
d$communeWaedenswil <- ifelse(d$commune == "Waedenswil", 1, 0)
d$communeWinterthur <- ifelse(d$commune == "Winterthur", 1, 0)
d$age31_50 <- ifelse(d$age_c3 == 1, 1, 0)
d$plus50 <- ifelse(d$age_c3 == 2, 1, 0)


#Estimate model
model1 <- glmmadmb(crossoff_pct ~nonswissnames + incumbent + relposition_rc + 
                      precumulated + sex_d + age31_50 + plus50 + CodeBerufK32 + CodeBerufK33 + num_cand + party_c32 + party_c33 +
                      communeAdliswil + communeBuelach + communeDietikon + communeWaedenswil + communeWinterthur + nonswissnames*relposition_rc +
                      (~1|list_id), data=d, family="Beta", admb.opts=admbControl(shess=FALSE,noinit=FALSE))

#-------------------- PLOT INTERACTION SWISS/NON-SWISS NAME & RELATIVE LIST POSITION

#Dataset for prediction
newdat <- expand.grid(
  crossoff_pct  = 0,
  relposition_rc  = seq(min(d$relposition_rc, rm=TRUE),
                        max(d$relposition_rc, rm=TRUE),
                        length.out=18),
  incumbent = 0, precumulated=0, sex_d=0, age31_50=0, plus50=0, party_c32=0, party_c33 = 0, num_cand=0,
  communeAdliswil = 0, communeBuelach=0, communeDietikon=0, communeWaedenswil=0, communeWinterthur=0, 
  CodeBerufK32=0, CodeBerufK33 = 0, nonswissnames = seq(0,1,1))

mm <- model.matrix(terms(model1), newdat)

#Calculate predictions

newdat$crossoff_norm <- predict(model1, 
                                newdat, re.form = NA, type="response")
pvar1 <- diag(mm %*% tcrossprod(vcov(model1), mm))


#Calculate lower/upper CI

cmult <- 1.96
newdat <- data.frame(
  newdat, plo = newdat$crossoff_norm - cmult*sqrt(pvar1), 
  phi = newdat$crossoff_norm + cmult*sqrt(pvar1))

#Plot of fixed effects uncertainty

newdat$Surname <- factor(newdat$nonswissnames,
                         levels = c(0,1),
                         labels = c("Swiss","Non-Swiss"))

marginaleffect <- ggplot(newdat, aes(relposition_rc, colour = Surname)) +
  geom_line(aes(y = crossoff_norm), size=1) +
  scale_color_manual(values = c("black", "gray60")) +
  geom_line(aes(y = plo,
                group=Surname), linetype =2, size=0.75) +
  geom_line(aes(y = phi,
                group=Surname), linetype =2, size=0.75) +
  xlab("Relative position on ballot") +
  ylab("Predicted values of negative preference votes") +
  labs(linetype='Non-Swiss surname') +
  scale_x_continuous(limits = c(-1, 0)) +
  scale_y_continuous(limits = c(-0.5, 1)) +
  theme_minimal() + 
  annotate("text", x = -0.89, y = -0.2, label = "Lowest position") +
  annotate("text", x = -0.12, y = -0.2, label = "Highest position")

#Layout
marginaleffect <- marginaleffect + theme(axis.text.y= element_text(size = 12, colour = "black")) + theme(axis.text.x= element_text(size = 12, colour = "black")) + 
  theme(axis.title = element_text(size = 12)) + theme(legend.title = element_text(size = 12)) + theme(legend.text = element_text(family="Times", size = 12)) + 
  theme(plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm"))
marginaleffect

tiff("Figure 2.tiff", width = 8.4, height = 6, units = 'in', res = 300)
plot(marginaleffect) # Make plot
dev.off()
